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Suppose we have a dataframe with n observations on m variables, i.e. m vectors in IR" . We want to 
transform the m variables in such a way that a particular aspect of the correlation matrix of the variables is 
maximized. An aspect is a real valued function/ defined on the space of correlation matrices. 

We formalize this using the framework discussed in De Leeuw (1988a; 1988b). Variable j defines a convex 
cone JCj in IR” . We suppose all vectors in the cone JCj are centered. Also define the sphere S of all vectors in 
IR” with sum of squares equal to one. For any choice of the Xj E JCj n S we can compute the correlation 
matrix R(X ) = {r^(X)} = xjx. ^. Our optimization problem is to maximize f(R(X)) over all Xj E JCj fi S. 

In our software we define the cones JCj as the sets of all vectors x such that x = Zb and (optionally) Ax > 0, 
where Z and A are gven matrices. In Mair and De Leeuw (2010) the r package aspect is presented. It 
deals with the special case where the variables are categorical and the Z are indicator matrices (a.k.a. 
dummies). On page 20 in Mair and De Leeuw (2010) it says: 

Considering each variable as categorical is not very efficient when having many 
categories, as typically in the numerical case. Therefore, in a future update we 
will use splines to make it more efficient. 


Well, this is that future update. In this note we develop r software that handles the much more general case 
where the Z are B-spline bases. The inequality constraints can be used to require that the spline 
transformations in Zb are monotone with the data. This covers the categorical case, because spline bases of 
degree zero are indicator matrices. It also covers polynomial transformations by using splines without interior 
knots. And finally our software also allows for nominal variables, by not imposing monotonicity restrictions. 

We now use a combination of tangential majorization and block relaxation (De Leeuw 2015a, De Leeuw 
(2015b)) to construct an algorithm maximizing/. Consider the case in which/ is convex on the convex set of 
correlation matrices. Then 


f(R(X)) >f(R(Y )) + tr G(Y)(R(X) - R(Y)), 

where G(Y ) := Df(R(Y)). In a majorization iteration we optimize tr G(Y)R(X). For this we use cyclic block 
relaxation, scaling one variable at a time while keeping the others fixed at their current values. Because 
diag(7?(X)) = diag(R(T)) we solve subproblem j in the cycle by optimizing xjtj over xj E JCj n S, where tj 
is the target 


h : = 2 SjAY)x f . 



Of course tj depends on Y and the Xf with € ^ j. We find the new Xj by minimizing (tj — Zjbj)'(tj — Zjbj ) 

A 

over AjZjbj > 0 and by projecting the optimal solution x-, = Zjbj on the sphere S. We assume throughout 
that Xj is non-zero, i.e. the target is not in the polar cone. 

In our r implementation we have incorporated some of the more obvious aspect choices that are convex 
functions of the correlation matrix: the sum of the dominant p eigenvalues (principal component analysis), the 
squared multiple correlation of one variable with the others (multiple regression), the sum of all squared 
multiple correlations of one variable with the others (image analysis), the negative logarithm of the 
determinant, the sum of pth powers of the correlations (with suitable p), the sum of pXh powers of the 
absolute values of the correlations. Thus our implementation covers both non-linear principal component 
analysis (Young, Takane, and Leeuw 1978, Winsberg and Ramsay (1983), Koyak (1987), Linting et al. (2007), 
De Leeuw (2014)) and multiple regression with optimal scaling (Kruskal 1965, Young, Leeuw, and Takane 
(1976), Winsberg and Ramsay (1980), Van Der Kooij, Meulman, and Heiser (2007)). And much more, obviously. 

Code 

aspect <- 

function (data, knots, degrees, ordinal, afunc, eps = le-6, itmax = 100, verbose = 

1, ...) { 

m <- ncol (data) 
n <- nrow (data) 
itel <- 1 

bd <- basesData (data, knots, degrees, ordinal) 
tdata <- matrix (0, n, m) 
for (j in l:m) { 

tdata[,j] <- bd$x[[j]] 

} 

tdata <- apply (tdata, 2, function (z) 
z - mean (z)) 

tdata <- apply (tdata, 2, function (z) 
z / sqrt (sum (z A 2))) 
corr <- crossprod (tdata) 
af <- afunc(corr, ...) 
fold <- af$f 
g <- af$g 
repeat { 

for (j in l:m) { 

target <- drop (tdata[,-j] %*% g[— j, j]) 
k <- bd$b[[j]] 
v <- bd$v[[j]] 

u <- drop (crossprod(k, target)) 
sO <- sum(target * tdata[,j]) 
if (ordinal[j]) { 
ns <- nnls(v,u) 
rr <- residuals(ns) 
tt <- drop(k %*% rr) 

} else 

tt <- drop (k %*% u) 
tt <- tt - mean (tt) 



sq <- sum(tt A 2) 
if (sq > le-15) { 

tt <- tt / sqrt (sum (tt A 2)) 
tdata[,j] <- tt 

} 

si <- sum(target * tdata[,j]) 
if (verbose > 1) 
cat ( 

"**** variable", formatC(j, width = 3, format = "d"), 
"Before", formate( 

sO, digits = 8, width = 12, format = "f" 

) r 

"After", formate( 

si, digits = 8, width = 12, format = "f" 

), 

" \ n" 

) 

corr <- crossprod (tdata) 
af <- afunc (corr, •••) 
fnew <- af$f 
g <- af$g 

} 

if (verbose > 0) 
cat ( 

"Iteration: ", formate (itel, width = 3, format = "d"), 
"fold: ", formate ( 

fold, digits = 8, width = 12, format = "f" 

) r 

"fnew: ", formate ( 

fnew, digits = 8, width = 12, format = "f" 

), 



if ((itel == itmax) || ((fnew - fold) < eps)) 

break 

itel <- itel + 1 
fold <- fnew 

} 

return (list ( 

tdata = tdata, f = fnew, r = corr, g = g 

)) 

} 

The program uses the nnis package (Mullen and Stokkum 2012) for the inequality restricted regressions, 
with the dual algorithm described by De Leeuw (2015c). The inequality constraints always require 
monotonicity of the transformed values Zb with the original data. 

The function aspect starts out by calling the subroutine basesData , which constructs and stores 
orthonormalized B-spline bases, initial transformation estimates, and other auxilary quantities in lists. The 
subroutine uses an ad-hoc B-spline implementation written in c (De Leeuw 2015d), and a similarly ad-hoc 



Gram-Schmidt in c (De Leeuw 2015e). For reproducibility purposes the c code is incorporated in this 
document. Note that we do not use monotone splines, such as the integrated B-splines of Winsberg and 
Ramsay (1980, -@winsberg_ramsay_83 (mailto:-@winsberg_ramsay_83)), we merely require that our splines 
are monotone at the data points. 

basesData <- function (data, knots, degrees, ordinal) { 
m <- ncol (data) 
n <- nrow (data) 
b <- as.list (rep (0, m)) 
v <- as.list (rep (0, m)) 
x <- as.list (rep (0, m)) 
for (j in l:m) { 

b[ [ j ] ] <- bsplineBasis(data[[j]], degrees[j], knots[[j]]) 
nb <- ncol(b[[j]]) 

X[[j]] <- drop (b[[j]] %*% (1 : nb)) 
b[[j]] <- apply (b[[j] ] , 2, function (z) 
z - mean (z)) 

b[[j]] <- (gsrc (b[[j]])$q)[,1 : (nb - 1), drop=FALSE] 
if (ordinal[j]) { 

v[ [ j ] ] <- matrix (0, nb - 1, n - 1) 
r <- rank (data[[j]], ties.method = "first") 
for (iinl : (n - 1)) { 

V[[ j]][, i] <- b[[j]][which(r == i),] - b[[j]][which(r == (i + 1)),] 

} 

} 

} 

return (list (b = b, v = v, x = x)) 

} 

Here are the implementations of various aspects. Of course the user can add her own aspect by writing a 
similar function. 

maxeig <- function (r, p) { 
e <- eigen (r) 
f <- sum (e$values[1:p]) 
g <- tcrossprod(e$vectors[,1:p]) 
return (list (f = f, g = g)) 

} 

maxcor <- function (r, p) { 
f <- sum (r A p) 
g <- P * (r * (p - 1)) 
return (list (f = f, g = g)) 

} 

maxabs <- function (r, p) { 
f <- sum (abs(r) A p) 
g <- P * (abs(r) " (p - 1))*sign(r) 
return (list (f = f, g = g)) 




} 


maxdet <- function (r) { 
f <- -log(det (r)) 
g <- -solve(r) 

return (list (f = f, g = g)) 

} 

maxsmc <- function (r, p) { 

beta <- solve (r[-p, -p], r[p, -p]) 

f <- sum (beta * r[p, -p]) 

h <- rep (1, nrow (r)) 

h[-p] <- -beta 

g <- -outer (h, h) 

return (list (f = f, g = g)) 

} 

maxsum <- function (r , p) { 
f <- sum (sqrt (r A 2 + p)) 
g <- r / sqrt (r A 2 + p) 
return (list (f = f, g = g)) 

} 

maximage <- function (r) { 
n <- nrow(r) 
f <- 0 

g <- matrix (0, n, n) 
for (p in l:n) { 

beta <- solve (r[-p,-p], r[p,-p]) 
f <- f + sum (beta * r[p,-p]) 
h <- rep (1, nrow (r)) 
h[-p] <- -beta 
g <- g - outer (h, h) 

} 

return (list (f = f, g = g)) 

} 

Example 1: Gibbs-Neumann Data 

Our first example are the Neumann data used by Willard Gibbs in his study of mixtures of gases with 
convertible components. These data have also been analyzed in Gifi (1990, 370-75), where there is more 
information about them. We choose the three knots at the quartiles. 

data(neumann, package = "homals") 

neumann_knots <- lapply (neumann, function (x) fivenum (x)[2:4]) 

neumann_degrees <- c(0,2,2) 

neumann_ordinal <-c(FALSE, TRUE , TRUE) 




First we maximize the multiple correlation of the last variable (density) with the first two (temperature and 
pressure). Temperature is scaled as nominal and categorical, pressure and density are ordinal with spline 
degree 2. 


hreg<-aspect(neumann,neumann_knots,neumann_degrees,neumann_ordinal , maxsmc,p=3) 


## 

Iteration: 

1 

fold: 

0.88154151 

f new: 

0.89539023 

## 

Iteration: 

2 

fold: 

0.89539023 

f new: 

0.89566693 

## 

Iteration: 

3 

fold: 

0.89566693 

f new: 

0.89566976 

## 

Iteration: 

4 

fold: 

0.89566976 

f new: 

0.89567005 


This gives the correlation matrix 


## 


[,1] 

[,2] 

[,3] 

## 

[i/i 

1.0000000 

0.3169096 

-0.8141950 

## 

[2,] 

0.3169096 

1.0000000 

0.1995548 

## 

[3,] 

-0.8141950 

0.1995548 

1.0000000 


and the optimal transformations (with vertical lines at the knots) 
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Next we compute a one-dimensional PCA by maximizing the largest eigenvalue (using the same knots and 
degrees). 


hpca<-aspect(neumann,neumann_knots,neumann_degrees,neumann_ordinal , maxeig,p=l) 


## 

Iteration: 

1 

fold: 

1.83957300 

f new: 

1.91034190 

## 

Iteration: 

2 

fold: 

1.91034190 

f new: 

1.91058869 

## 

Iteration: 

3 

fold: 

1.91058869 

f new: 

1.91059236 

## 

Iteration: 

4 

fold: 

1.91059236 

f new: 

1.91059268 


The correlation matrix now is 


## 


[,1] 

[,2] 

[,3] 

## 

[i# i 

1.0000000 

0.400810296 

-0.819241802 

## 

[2, ] 

0.4008103 

1.000000000 

0.003661109 

## 

[3, ] 

-0.8192418 

0.003661109 

1.000000000 


and the optimal transformations are 
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Next, we maximize what is perhaps the simplest aspect: the sum of the correlations. 


hsum<-aspect(neumann,neumann_knots,neumann_degrees,neumann_ordinal,maxcor,p=l) 


## 

Iteration: 

1 

fold: 

2.40147297 

f new: 

4.20649363 

## 

Iteration: 

2 

fold: 

4.20649363 

f new: 

4.20663743 

## 

Iteration: 

3 

fold: 

4.20663743 

f new: 

4.20669120 

## 

Iteration: 

4 

fold: 

4.20669120 

f new: 

4.20671222 

## 

Iteration: 

5 

fold: 

4.20671222 

f new: 

4.20672045 

## 

Iteration: 

6 

fold: 

4.20672045 

f new: 

4.20672367 

## 

Iteration: 

7 

fold: 

4.20672367 

f new: 

4.20672493 

## 

Iteration: 

8 

fold: 

4.20672493 

f new: 

4.20672543 


This gives the correlation matrix 


hsum$r 


## 


[,1] 

[,2] 

[ / 3 ] 

## 

[i/i 

1.0000000 

-0.3596036 

0.8030536 

## 

[2,] 

-0.3596036 

1.0000000 

0.1599128 

## 

[3,] 

0.8030536 

0.1599128 

1.0000000 










The sum of the absolute values of the correlations. 


habs<-aspect(neumann,neumann_knots,neumann_degrees,neumann_ordinal,maxabs,p=l) 


## 

Iteration: 

1 

fold: 

5.64234387 

f new: 

5.66909477 

## 

Iteration: 

2 

fold: 

5.66909477 

f new: 

5.66992570 

## 

Iteration: 

3 

fold: 

5.66992570 

f new: 

5.66996738 

## 

Iteration: 

4 

fold: 

5.66996738 

f new: 

5.66997382 

## 

Iteration: 

5 

fold: 

5.66997382 

f new: 

5.66997485 

## 

Iteration: 

6 

fold: 

5.66997485 

f new: 

5.66997501 


In this case the correlation matrix after the last iteration is 


## 


[,1] 

[,2] 

[,3] 

## 

[i# i 

1.0000000 

0.3209786 

-0.8013430 

## 

[2, ] 

0.3209786 

1.0000000 

0.2126659 

## 

[3, ] 

-0.8013430 

0.2126659 

1.0000000 


and the optimal transformations are 
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Example 2: Bodyfat Data 

Age, weight, height, and 10 body circumference measurements are recorded for 252 men. Each man’s 
percentage of body fat was accurately estimated by an underwater weighing technique. 


data(fat, package = "faraway") 
bodyfat <- fat 

bodyfat_knots <- lapply (bodyfat, function (x) fivenum (x)[2:4]) 

bodyfat_degrees <- rep(2, 18) 

bodyfat_ordinal <- rep (TRUE, 18) 

bodyfat_ordinal[4]<-FALSE 

names(bodyfat) 


## 

[1] 

"brozek" 

"siri" 

"density" 

"age" 

"weight" 

"height" 

"adipos" 

## 

[8] 

"free" 

"neck" 

"chest" 

"abdom" 

"hip" 

"thigh" 

"knee" 

## 

[15] 

"ankle" 

"biceps" 

"forearm" 

"wrist" 





The first two variables are bodyfat indices, computed with different formulas. We see how well we can predict 
the first one using the multiple correlation aspect. All variables are transformed monotonically, except age. 


bsmc<-aspect(bodyfat,bodyfat_knots,bodyfat_degrees,bodyfat_ordinal,maxsmc,p=l) 









## 

Iteration: 

1 

fold 

0.99937103 

f new: 

0.99971962 

## 

Iteration: 

2 

fold 

0.99971962 

f new: 

0.99974700 

## 

Iteration: 

3 

fold 

0.99974700 

f new: 

0.99975850 

## 

Iteration: 

4 

fold 

0.99975850 

f new: 

0.99976400 

## 

Iteration: 

5 

fold 

0.99976400 

f new: 

0.99976708 

## 

Iteration: 

6 

fold 

0.99976708 

f new: 

0.99976926 

## 

Iteration: 

7 

fold 

0.99976926 

f new: 

0.99977107 

## 

Iteration: 

8 

fold 

0.99977107 

f new: 

0.99977264 

## 

Iteration: 

9 

fold 

0.99977264 

f new: 

0.99977406 

## 

Iteration: 

10 

fold 

0.99977406 

f new: 

0.99977537 

## 

Iteration: 

11 

fold 

0.99977537 

f new: 

0.99977659 

## 

Iteration: 

12 

fold 

0.99977659 

f new: 

0.99977775 

## 

Iteration: 

13 

fold 

0.99977775 

f new: 

0.99977885 

## 

Iteration: 

14 

fold 

0.99977885 

f new: 

0.99977990 

## 

Iteration: 

15 

fold 

0.99977990 

f new: 

0.99978090 

## 

Iteration: 

16 

fold 

0.99978090 

f new: 

0.99978184 


We plot the optimal transformation for the brozek index and for age. 
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We can use nonlinear principal component analysis to look more closely at the structure of the variables 
defining the indices. Optimize the sum of the two largest eigenvalues of the correlation matix of the last 16 
variables. 


beig<-aspect(bodyfat[,3:18], bodyfat_knots[3:18], bodyfat_degrees[3:18], bodyfat_ordi 
nal[3:18], maxeig, p = 2) 







## 

Iteration: 

1 

fold: 

11.93445666 

f new: 

12.23316284 

## 

Iteration: 

2 

fold: 

12.23316284 

f new: 

12.24922704 

## 

Iteration: 

3 

fold: 

12.24922704 

f new: 

12.25509701 

## 

Iteration: 

4 

fold: 

12.25509701 

f new: 

12.25870256 

## 

Iteration: 

5 

fold: 

12.25870256 

f new: 

12.26121744 

## 

Iteration: 

6 

fold: 

12.26121744 

f new: 

12.26300909 

## 

Iteration: 

7 

fold: 

12.26300909 

f new: 

12.26429152 

## 

Iteration: 

8 

fold: 

12.26429152 

f new: 

12.26521099 

## 

Iteration: 

9 

fold: 

12.26521099 

f new: 

12.26587082 

## 

Iteration: 

10 

fold: 

12.26587082 

f new: 

12.26634466 

## 

Iteration: 

11 

fold: 

12.26634466 

f new: 

12.26668514 

## 

Iteration: 

12 

fold: 

12.26668514 

f new: 

12.26692993 

## 

Iteration: 

13 

fold: 

12.26692993 

f new: 

12.26710602 

## 

Iteration: 

14 

fold: 

12.26710602 

f new: 

12.26723274 

## 

Iteration: 

15 

fold: 

12.26723274 

f new: 

12.26732397 

## 

Iteration: 

16 

fold: 

12.26732397 

f new: 

12.26738968 

## 

Iteration: 

17 

fold: 

12.26738968 

f new: 

12.26743702 

## 

Iteration: 

18 

fold: 

12.26743702 

f new: 

12.26747113 

## 

Iteration: 

19 

fold: 

12.26747113 

f new: 

12.26749573 

## 

Iteration: 

20 

fold: 

12.26749573 

f new: 

12.26751346 

## 

Iteration: 

21 

fold: 

12.26751346 

f new: 

12.26752624 

## 

Iteration: 

22 

fold: 

12.26752624 

f new: 

12.26753546 

## 

Iteration: 

23 

fold: 

12.26753546 

f new: 

12.26754212 

## 

Iteration: 

24 

fold: 

12.26754212 

f new: 

12.26754692 

## 

Iteration: 

25 

fold: 

12.26754692 

f new: 

12.26755038 

## 

Iteration: 

26 

fold: 

12.26755038 

f new: 

12.26755288 

## 

Iteration: 

27 

fold: 

12.26755288 

f new: 

12.26755468 

## 

Iteration: 

28 

fold: 

12.26755468 

f new: 

12.26755598 

## 

Iteration: 

29 

fold: 

12.26755598 

f new: 

12.26755692 


The eigenvalues of the correlation matrix, in percentages, with a plot of their cumulative values, are as follows. 


## 

[i] 

0.6369047237 

0.1298175838 

0.0593329100 

0.0399460439 

0.0308797444 

## 

[6] 

0.0285711710 

0.0181510533 

0.0148222809 

0.0128746126 

0.0110294094 

## 

[11] 

0.0076337401 

0.0042483034 

0.0026208110 

0.0021778757 

0.0006752046 

## 

[16] 

0.0003145321 
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The scatterplot of the first two principal components is as follows. 
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Example 3: Air Pollution 

The Sokal/Rolph air pollution data have previously been used by Wang and Murphy (2004) to illustrate ACE 
(Breiman and Friedman 1985). The outcome variable is sulpher dioxide concentration in air mgs. per cubic 
metre in 41 cities in the USA. The predictors are 

• average annual temperature in degrees F 

• number of manufacturers employing >20 workers 

• population size in thousands 

• average annual wind speed in miles per hour 

• average annual rainfall in inches 

• average number of days rainfall per year 


data(usair, package = "gamlss.data") 

usair_knots <- lapply (usair, function (x) fivenum (x)[2:4]) 
usair_degrees <- c (3,3,3,3,3,3,3) 

USair_ordinal <- C ( TRUE , TRUE , TRUE , TRUE , FALSE , TRUE , FALSE ) 


airsmc<-aspect(usair,usair_knots,usair_degrees,usair_ordinal,verbose=0, itmax = 1000, 
maxsmc,p=l) 







The optimal SMC is 0.9482315 (for ACE this is 0.9469), the correlation matrix is 


## 


[,1] 

[,2] 

[,3] 

[,4] 

[,5] 

[,6] 

## 

[1#] 

1.0000000 

-0.34401533 

0.85100005 

0.34248610 

-0.4040813 

0.04775350 

## 

[2,] 

-0.3440153 

1.00000000 

-0.16310443 

0.06765795 

0.1208380 

0.34046001 

## 

[3,] 

0.8510000 

-0.16310443 

1.00000000 

0.69790740 

-0.3560280 

-0.03922927 

## 

[4,] 

0.3424861 

0.06765795 

0.69790740 

1.00000000 

-0.3947447 

-0.05441145 

## 

[5,] 

-0.4040813 

0.12083797 

-0.35602805 

-0.39474473 

1.0000000 

0.07234770 

## 

[6,] 

0.0477535 

0.34046001 

-0.03922927 

-0.05441145 

0.0723477 

1.00000000 

## 

[7,] 

-0.1707223 

-0.24678998 

-0.09406506 

-0.13570731 

-0.1345694 

-0.05539428 

## 


[,7] 






## 

[1#] 

-0.17072228 






## 

[2,] 

-0.24678998 






## 

[3,] 

-0.09406506 






## 

[4,] 

-0.13570731 






## 

[5,] 

-0.13456940 






## 

[6,] 

-0.05539428 






## 

[7,] 

1.00000000 






and the regression coefficients are 

## 

[1] 

-0.2058751 

1.0722221 -0 

.5066085 -0. 

2361394 0. 

1375961 -0. 

2135771 


The seven transformation plots are as follows. 
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We also use this example to show that an empty vector of interior knots leads to a polynomial transformation. 
In this case we transform the dependent variable linearly (which basically means we do not transform it at all). 


usair_knots[[1]]<-numeric(0) 
usair_degrees[1]<-l 

airlin<-aspect(usair,usair_knots,usair_degrees,usair_ordinal,verbose=0,maxsmc, p=l) 


The optimal SMC is now 0.904187, the correlation matrix is 



























## 


[,1] 

[,2] 

[,3] 

[,4] 

[,5] 

## 

[1, ] 

1.000000000 

-0.4829229 

0.60609579 

0.28419105 

-0.36590731 

## 

[2, ] 

-0.482922926 

1.0000000 

-0.19383688 

0.01929680 

0.15465911 

## 

[3, ] 

0.606095786 

-0.1938369 

1.00000000 

0.86739171 

-0.35180667 

## 

[4, ] 

0.284191054 

0.0192968 

0.86739171 

1.00000000 

-0.37329140 

## 

[5, ] 

-0.365907315 

0.1546591 

-0.35180667 

-0.37329140 

1.00000000 

## 

[6, ] 

0.175496105 

0.2817608 

-0.02235636 

-0.03976988 

-0.01936234 

## 

[7, ] 

0.002801162 

-0.3629289 

0.05071848 

-0.03275958 

-0.22797737 

## 


[,6] 

[ r 7] 




## 

[1, ] 

0.17549611 

0.002801162 




## 

[2, ] 

0.28176080 

-0.362928895 




## 

[3, ] 

-0.02235636 

0.050718480 




## 

[4, ] 

-0.03976988 

-0.032759583 




## 

[5, ] 

-0.01936234 

-0.227977370 




## 

[6, ] 

1.00000000 

0.359290635 




## 

[7, ] 

0.35929063 

1.000000000 





and the regression coefficients are 


## [1] -0.5412443 1.1010932 -0.7656825 -0.2898311 0.5039996 -0.5217185 


The seven transformation plots are as follows. 
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Example 4: Moral Integration of American Cities 

The angell data frame has 43 rows and 4 columns. The observations are 43 U. S. cities around 1950. These 
data were previously analyzed with ACE by De Veaux (1987). The variables are 

• Moral Integration: Composite of crime rate and welfare expenditures. 

• Ethnic Heterogenity: From percentages of nonwhite and foreign-born white residents. 

• Geographic Mobility: From percentages of residents moving into and out of the city. 

• A factor with levels: Northeast; Midwest; Southeast; West. 


data(Angell, package = "car") 
angell <- Angell[,1:3] 

angell_knots <- lapply (angell, function (x) fivenum (x)[2:4]) 
angell_degrees <- c(2,2,2) 
angell_ordinal <- c(TRUE, TRUE, TRUE) 


hsmc<-aspect(angell,angell_knots,angell_degrees,angell_ordinal,maxsmc,p=l) 




























## 

Iteration: 

1 

fold: 

0.67555675 

f new: 

0.73923854 

## 

Iteration: 

2 

fold: 

0.73923854 

f new: 

0.74680608 

## 

Iteration: 

3 

fold: 

0.74680608 

f new: 

0.74931438 

## 

Iteration: 

4 

fold: 

0.74931438 

f new: 

0.75005255 

## 

Iteration: 

5 

fold: 

0.75005255 

f new: 

0.75026068 

## 

Iteration: 

6 

fold: 

0.75026068 

f new: 

0.75031127 

## 

Iteration: 

7 

fold: 

0.75031127 

f new: 

0.75032348 

## 

Iteration: 

8 

fold: 

0.75032348 

f new: 

0.75032642 

## 

Iteration: 

9 

fold: 

0.75032642 

f new: 

0.75032713 


Correlations between transformed variables are as follows. 


## 


[,1] 

[,2] 

[,3] 

## 

[i# i 

1.0000000 

-0.53934872 

-0.64048623 

## 

[2, ] 

-0.5393487 

1.00000000 

-0.06643057 

## 

[3, ] 

-0.6404862 

-0.06643057 

1.00000000 


The transformation plots are as follows. 
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Example 5: Thirteen Personality Scales 

We illustrate a more complicated aspect using a small data set from the psych package (Revelle 2015) of five 
scales from the Eysenck Personality Inventory, five from a Big 5 inventory, a Beck Depression Inventory, and 
State and Trait Anxiety measures. 


data(epi.bfi, package = "psych") 
epi <- epi.bfi 

epi_knots <- lapply (epi, function (x) fivenum (x)[2:4]) 
epi_degrees <- rep (2, 13) 
epi_ordinal <- rep (TRUE, 13) 


The aspect we use is the maximum log-likelihood of a factor analysis model. Because this aspect is the 
maximum of a family of functions linear in r it is indeed convex. Obviously the same result applies to any 
structural equation model combined with a multinormal log-likelihood. Using the negative logarithm of the 
determinant as an aspect corresponds with the special case of a saturated model. 











maxfac <- function (r, p) { 

fa <- factanal (NULL, p, covmat = r, rotation = "none") 
s <- tcrossprod (fa$loadings) + diag (fa$unique) 
g <- - solve (s) 
f <- sum (g * r) - log(det (s)) 
return (list (f = f, g = g)) 

> 


We fit a two-factor model, without worrying if two factors are appropriate for these data. 


efac<-aspect(epi,epi_knots,epi_degrees,epi_ordinal,maxfac,p=2) 


## 

Iteration: 

1 

fold 

-7.47534413 

f new: 

-7.08355786 

## 

Iteration: 

2 

fold 

-7.08355786 

f new: 

-7.05118236 

## 

Iteration: 

3 

fold 

-7.05118236 

f new: 

-7.03813254 

## 

Iteration: 

4 

fold 

-7.03813254 

f new: 

-7.03277037 

## 

Iteration: 

5 

fold 

-7.03277037 

f new: 

-7.03050114 

## 

Iteration: 

6 

fold 

-7.03050114 

f new: 

-7.02952915 

## 

Iteration: 

7 

fold 

-7.02952915 

f new: 

-7.02911088 

## 

Iteration: 

8 

fold 

-7.02911088 

f new: 

-7.02893056 

## 

Iteration: 

9 

fold 

-7.02893056 

f new: 

-7.02885277 

## 

Iteration: 

10 

fold 

-7.02885277 

f new: 

-7.02881920 

## 

Iteration: 

11 

fold 

-7.02881920 

f new: 

-7.02880472 

## 

Iteration: 

12 

fold 

-7.02880472 

f new: 

-7.02879847 

## 

Iteration: 

13 

fold 

-7.02879847 

f new: 

-7.02879577 

## 

Iteration: 

14 

fold 

-7.02879577 

f new: 

-7.02879461 

## 

Iteration: 

15 

fold 

-7.02879461 

f new: 

-7.02879411 


The transformation plots are as follows. 
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Here are the factor loadings, for those so inclined 


















o 


LO 

o 


CNJ 

o q 

o O 

.03 


LO 

O 



factor 1 


1.5 


And here are the uniquenesses. Variable 1 (epiE) seems to be heading for Heywood. 


## 

[1] 

0.0050000 

0.2527709 

0.3638904 

0.8268849 

0.3507721 

0.8915067 0.8864692 

## 

[8] 

0.6280157 

0.5196517 

0.9219710 

0.4660321 

0.1194840 

0.5526974 
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